Stephens lab data accessed from http://www.gtexportal.org/static/datasets/gtex_analysis_pilot_v3/multi_tissue_eqtls/Multi_tissue_eQTL_GTEx_Pilot_Phase_datasets.tar on 20150722.
From README: “We are using the eQTL posterior probabilities from the UC Multi-tissue eQTL method (doi:10.1371/journal.pgen.1003486) for each of the 9 tissues analyzed in the pilot phase (Adipose_subcutaneous, Artery_Tibal, Whole_Blood, Heart_Left_Ventricle, Lung, Muscle_Skeletal, Nerve_Tibial, Skin_Lower_Leg_Sun_Exposed, Thyroid) in the file res_final_uc_com_genes_com_snps.txt.gz. These values may be interpreted as Pr(SNP is eQTL in tissue s | data). 9875 eGenes are presented, with the”top" (most significant) SNP in each gene used."
library(plyr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(GGally)
"%&%" = function(a,b) paste(a,b,sep="")
mt.dir <- "~/GitHub/GenArch/GenArchPaper/UC_Stephens_Multitissue_Comp/"
#tarbell mt.dir <- "/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/Multi_tissue_eQTL_GTEx_Pilot_Phase_datasets/"
h2.dir <- "~/GitHub/GenArch/GenArchPaper/BSLMM/bslmm_gtex_results/"
#tarbell h2.dir <- "/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/"
bslmm.dir <- "~/GitHub/GenArch/GenArchPaper/BSLMM/bslmm_gtex_results/"
#tarbell bslmm.dir <- "/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/BSLMM_exp/"
mt <- read.table(mt.dir %&% "res_final_uc_com_genes_com_snps.txt.gz",header=TRUE)
head(mt)
## gene snp Adipose Artery Blood Heart Lung Muscle
## 1 ENSG00000137288.4 rs11759971 0.6596 0.6547 0.5383 0.6264 0.6759 0.5699
## 2 ENSG00000168005.4 rs36077346 0.7835 0.7817 0.7632 0.7743 0.7928 0.5103
## 3 ENSG00000255893.1 rs657249 1.0000 1.0000 0.3647 1.0000 1.0000 1.0000
## 4 ENSG00000103227.12 rs2382946 0.9812 0.9985 0.9967 1.0000 0.9998 1.0000
## 5 ENSG00000188385.7 rs9733324 0.9939 0.9935 0.6668 0.9839 0.9949 0.9801
## 6 ENSG00000166257.4 rs1148088 0.4582 0.4515 0.3038 0.3608 0.5185 0.2841
## Nerve Skin Thyroid
## 1 0.6750 0.7638 0.6775
## 2 0.7922 0.7837 0.7932
## 3 1.0000 1.0000 1.0000
## 4 0.9988 0.9104 1.0000
## 5 0.9945 0.9926 0.9950
## 6 0.4813 0.4602 0.8391
dim(mt)
## [1] 9875 11
#remove version number in order to compare ensembl IDs
a <- substr(mt$gene,1,15)
mt <- mutate(mt,gene=a)
h2.ts <- read.table(bslmm.dir %&% "GTEx_Tissue-Specific_local_h2.txt",header=TRUE)
#remove version number in order to compare ensembl IDs
a <- substr(h2.ts$ensid,1,15)
h2.ts <- mutate(h2.ts,gene=a)
h2.tw <- read.table(bslmm.dir %&% "GTEx_Tissue-Wide_local_h2.txt",header=TRUE)
#remove version number in order to compare ensembl IDs
a <- substr(h2.tw$ensid,1,15)
h2.tw <- mutate(h2.tw,gene=a)
##Get BSLMM results
pve.ts <- read.table(bslmm.dir %&% "GTEx_Tissue-Specific_local_PVE_by_BSLMM.txt",header=TRUE)
#remove version number in order to compare ensembl IDs
pve.ts <- mutate(pve.ts,gene=substr(pve.ts$gene,1,15))
pve.tw <- read.table(bslmm.dir %&% "GTEx_Tissue-Wide_local_PVE_by_BSLMM.txt",header=TRUE)
#remove version number in order to compare ensembl IDs
pve.tw <- mutate(pve.tw,gene=substr(pve.tw$gene,1,15))
entropy <- function(pr){
pi <- pr/sum(pr)
pi[pi==0] <- 1e-06 ##allows log transformation
-1*sum(pi*log(pi))
}
mtPr <- mt[,3:11]
mtS <- apply(mtPr,1,entropy) ##calc entropy for each row
hist(mtS)
summary(mtS)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002125 2.149000 2.183000 2.127000 2.196000 2.197000
geneS <- select(mt,gene) %>% mutate(entropy=mtS) ##put gene and entropy together
tislist <- c('CrossTissue','AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
h2.tis <- h2.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(geneS,h2.tis,by='gene')
colnames(newdata)[3] <- "h2"
p <- ggplot(newdata,aes(x=h2,y=entropy))+geom_point(alpha=0.4)+stat_smooth()+xlab(tis %&% " h2")
cors <- data.frame(cor=signif(cor(newdata$h2, newdata$entropy), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=4)
pvals <- data.frame(cor=signif(cor.test(newdata$h2, newdata$entropy)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=4)
print(p2)
}
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
h2.tis <- h2.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(geneS,h2.tis,by='gene')
colnames(newdata)[3] <- "h2"
p <- ggplot(newdata,aes(x=h2,y=entropy))+geom_point(alpha=0.4)+stat_smooth()+
xlab(tis %&% " tissue-wide h2")
cors <- data.frame(cor=signif(cor(newdata$h2, newdata$entropy), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=4)
pvals <- data.frame(cor=signif(cor.test(newdata$h2, newdata$entropy)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=4)
print(p2)
}
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
tislist <- c('CrossTissue','AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
pve.tis <- pve.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(geneS,pve.tis,by='gene')
colnames(newdata)[3] <- "pve"
p <- ggplot(newdata,aes(x=pve,y=entropy))+geom_point(alpha=0.4)+stat_smooth()+
xlab(tis %&% " PVE")
cors <- data.frame(cor=signif(cor(newdata$pve, newdata$entropy, use="p"), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=4)
pvals <- data.frame(cor=signif(cor.test(newdata$pve, newdata$entropy)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=4)
print(p2)
}
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
pve.tis <- pve.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(geneS,pve.tis,by='gene')
colnames(newdata)[3] <- "pve"
p <- ggplot(newdata,aes(x=pve,y=entropy))+geom_point(alpha=0.4)+stat_smooth()+
xlab(tis %&% " tissue-wide PVE")
cors <- data.frame(cor=signif(cor(newdata$pve, newdata$entropy, use="p"), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=4)
pvals <- data.frame(cor=signif(cor.test(newdata$pve, newdata$entropy)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=4)
print(p2)
}
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
data <- inner_join(mt,h2.ts,by='gene')
dim(data)
## [1] 7173 22
##only half overlap, why? gcta didn't converge? multi-tissue model didn't converge?
#test one tissue h2 at a time
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
h2.tis <- h2.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(mt,h2.tis,by='gene')
gdata <- gather(newdata,"tissue","Pr",3:11)
colnames(gdata)[3] <- "h2"
p <- ggplot(gdata,aes(x=h2,y=Pr))+facet_wrap(~tissue)+geom_point(alpha=0.4)+stat_smooth()+
coord_cartesian(ylim=c(-0.05,1.2))+xlab(tis %&% " h2")
cors <- ddply(gdata, .(tissue), summarise, cor = signif(cor(h2, Pr), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=3)
pvals <- ddply(gdata, .(tissue), summarise, cor = signif(cor.test(h2, Pr)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=3)
print(p2)
}
tis <- 'CrossTissue'
h2.tis <- h2.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(mt,h2.tis,by='gene')
gdata <- gather(newdata,"tissue","Pr",3:11)
colnames(gdata)[3] <- "h2"
p <- ggplot(gdata,aes(x=h2,y=Pr))+facet_wrap(~tissue)+geom_point(alpha=0.4)+stat_smooth()+
coord_cartesian(ylim=c(-0.05,1.2))+xlab(tis %&% " h2")
cors <- ddply(gdata, .(tissue), summarise, cor = signif(cor(h2, Pr), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=3)
pvals <- ddply(gdata, .(tissue), summarise, cor = signif(cor.test(h2, Pr)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=3)
print(p2)
#test one tissue h2 at a time
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
h2.tis <- h2.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(mt,h2.tis,by='gene')
gdata <- gather(newdata,"tissue","Pr",3:11)
colnames(gdata)[3] <- "h2"
p <- ggplot(gdata,aes(x=h2,y=Pr))+facet_wrap(~tissue)+geom_point(alpha=0.4)+stat_smooth()+
coord_cartesian(ylim=c(-0.05,1.2))+xlab(tis %&% " tissue-wide h2")
cors <- ddply(gdata, .(tissue), summarise, cor = signif(cor(h2, Pr), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=3)
pvals <- ddply(gdata, .(tissue), summarise, cor = signif(cor.test(h2, Pr)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=3)
print(p2)
}
data <- inner_join(mt,pve.ts,by='gene')
dim(data)
## [1] 7173 21
##only half overlap, why? same overlap as with h2? max overlap b/t gencode versions?
#test one tissue pve at a time
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
pve.tis <- pve.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(mt,pve.tis,by='gene')
gdata <- gather(newdata,"tissue","Pr",3:11)
colnames(gdata)[3] <- "pve"
p <- ggplot(gdata,aes(x=pve,y=Pr))+facet_wrap(~tissue)+geom_point(alpha=0.4)+stat_smooth()+
coord_cartesian(ylim=c(-0.05,1.2))+xlab(tis %&% " PVE")
cors <- ddply(gdata, .(tissue), summarise, cor = signif(cor(pve, Pr, use='p'), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=3)
pvals <- ddply(gdata, .(tissue), summarise, cor = signif(cor.test(pve, Pr)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=3)
print(p2)
}
tis <- 'CrossTissue'
pve.tis <- pve.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(mt,pve.tis,by='gene')
gdata <- gather(newdata,"tissue","Pr",3:11)
colnames(gdata)[3] <- "pve"
p <- ggplot(gdata,aes(x=pve,y=Pr))+facet_wrap(~tissue)+geom_point(alpha=0.4)+stat_smooth()+
coord_cartesian(ylim=c(-0.05,1.2))+xlab(tis %&% " PVE")
cors <- ddply(gdata, .(tissue), summarise, cor = signif(cor(pve, Pr, use="p"), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=3)
pvals <- ddply(gdata, .(tissue), summarise, cor = signif(cor.test(pve, Pr)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=3)
print(p2)
#test one tissue pve at a time
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
pve.tis <- pve.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(mt,pve.tis,by='gene')
gdata <- gather(newdata,"tissue","Pr",3:11)
colnames(gdata)[3] <- "pve"
p <- ggplot(gdata,aes(x=pve,y=Pr))+facet_wrap(~tissue)+geom_point(alpha=0.4)+stat_smooth()+
coord_cartesian(ylim=c(-0.05,1.2))+xlab(tis %&% " tissue-wide PVE")
cors <- ddply(gdata, .(tissue), summarise, cor = signif(cor(pve, Pr, use="p"), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=3)
pvals <- ddply(gdata, .(tissue), summarise, cor = signif(cor.test(pve, Pr)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=3)
print(p2)
}
tislist <- c('CrossTissue','AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
pve.tis <- pve.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
h2.tis <- h2.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(h2.tis,pve.tis,by='gene')
colnames(newdata)[2:3] <- c("h2","pve")
p <- ggplot(newdata,aes(x=h2,y=pve))+geom_point(alpha=0.4)+stat_smooth()+xlab(tis %&% " h2")+ylab(tis %&% " PVE")
cors <- data.frame(cor=signif(cor(newdata$pve, newdata$h2, use="p"), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.1, size=3)
pvals <- data.frame(cor = signif(cor.test(newdata$pve, newdata$h2)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.1, vjust=1.1, size=3)
print(p2)
}
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
pve.tis <- pve.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
h2.tis <- h2.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(h2.tis,pve.tis,by='gene')
colnames(newdata)[2:3] <- c("h2","pve")
p <- ggplot(newdata,aes(x=h2,y=pve))+geom_point(alpha=0.4)+stat_smooth()+xlab(tis %&% " tissue-wide h2")+ylab(tis %&% " tissue-wide PVE")
cors <- data.frame(cor=signif(cor(newdata$pve, newdata$h2, use="p"), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.1, size=3)
pvals <- data.frame(cor = signif(cor.test(newdata$pve, newdata$h2)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.1, vjust=1.1, size=3)
print(p2)
}
meanmt <- apply(mt[,3:11],1,mean)
sdmt <- apply(mt[,3:11],1,sd)
cv <- sdmt/meanmt
newmt<-mutate(mt,CV=cv)
tislist <- c('CrossTissue','AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
h2.tis <- h2.ts %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(newmt,h2.tis,by='gene')
colnames(newdata)[13] <- "h2"
p <- ggplot(newdata,aes(x=h2,y=CV))+geom_point(alpha=0.4)+stat_smooth()+
xlab(tis %&% " h2")+ylab("Coef of variation (sd/mean) of Pr")
cors <- data.frame(cor=signif(cor(newdata$h2, newdata$CV), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=4)
pvals <- data.frame(cor=signif(cor.test(newdata$h2, newdata$CV)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=4)
print(p2)
}
meanmt <- apply(mt[,3:11],1,mean)
sdmt <- apply(mt[,3:11],1,sd)
cv <- sdmt/meanmt
newmt<-mutate(mt,CV=cv)
tislist <- c('AdiposeSubcutaneous','ArteryTibial','HeartLeftVentricle','Lung','MuscleSkeletal','NerveTibial','SkinSunExposedLowerleg','Thyroid','WholeBlood')
for(tis in tislist){
h2.tis <- h2.tw %>% select(gene,one_of(tis)) ##one_of allows character vector
newdata <- inner_join(newmt,h2.tis,by='gene')
colnames(newdata)[13] <- "h2"
p <- ggplot(newdata,aes(x=h2,y=CV))+geom_point(alpha=0.4)+stat_smooth()+
xlab(tis %&% " tissue-wide h2")+ylab("Coef of variation (sd/mean) of Pr")
cors <- data.frame(cor=signif(cor(newdata$h2, newdata$CV), 2))
p1 <- p + geom_text(data=cors, aes(label=paste("R = ", cor, sep="")), x=-Inf, y=Inf, hjust=-0.2, vjust=1.2, size=4)
pvals <- data.frame(cor=signif(cor.test(newdata$h2, newdata$CV)$p.value, 2))
p2 <- p1 + geom_text(data=pvals, aes(label=paste("p = ", cor, sep="")), x=Inf, y=Inf, hjust=1.2, vjust=1.2, size=4)
print(p2)
}